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We perform Bayesian analysis of gravitational-wave signals from non-spinning, intermediate-mass black- 
hole binaries (IMBHBs) with observed total mass, Mobs, from SOM© to SOOM© and mass ratio 1-4 using 
advanced LIGO and Virgo detectors. We employ inspiral-merger-ringdown waveform models based on the 
effective-one-body formalism and include subleading modes of radiation beyond the leading (2, 2) mode. The 
presence of subleading modes increases signal power for inclined binaries and allows for improved accuracy 
and precision in measurements of the masses as well as breaking of degeneracies in distance, orientation and 
polarization. For low total masses. Mobs ^ 50M©, for which the inspiral signal dominates, the observed chirp 
mass Mlobs = Mobs (j] being the symmetric mass ratio) is better measured. In contrast, as increasing 
power comes from merger and ringdown, we find that the total mass Mobs has better relative precision than 
Mohs- Indeed, at high Mobs (> 300M©), the signal resembles a burst and the measurement thus extracts the 
dominant frequency of the signal that depends on Mobs- Depending on the binary’s inclination, at signal-to- 
noise ratio (SNR) of 12, uncertainties in Mobs can be as large as ~ 20-25% while uncertainties in Afobs are 
~ 50-60% in binaries with unequal masses (those numbers become ~ 17% versus ~ 22% in more symmetric 
mass-ratio binaries). Although large, those uncertainties in Mobs will establish the existence of IMBHs. We 
find that effective-one-body waveforms with subleading modes are essential to confirm a signal’s presence in 
the data, with calculated Bayesian evidences yielding a false alarm probability below 10“^ for SNR > 9 
in Gaussian noise. Our results show that gravitational-wave observations can offer a unique tool to observe 
and understand the formation, evolution and demographics of IMBHs, which are difficult to observe in the 
electromagnetic window. 

PACS numbers: 04.30.Tv, 04.70.Bw, 02.70.Uu, 04.80.Nn 
Keywords: gravitational waves, black holes, Bayesian inference 


I. INTRODUCTION 

Advanced interferometric gravitational-wave (GW) detec¬ 
tors LIGO and Virgo will be turned on in late 2015 (2016 for 
Virgo) and are expected to reach design sensitivity by 2019 
[1-4]. At design sensitivity these detectors will operate in 
the frequency range from 10 Hz to 1 kHz, with an almost flat 
sensitivity from 40 Hz to 1 kHz. In Fig. 1, we show for the ad¬ 
vanced LIGO-Virgo network, the distance reach ^ as a function 
of observed total mass for the full inspiral-merger-ringdown 
signal of binaries consisting of nonspinning black holes (BHs) 
and include several multipole modes beyond the dominant 
(^ = 2, m = 2) mode, as well as higher order post-Newtonian 
(PN) corrections (see Sec. II for details) [5]. For nonspin¬ 
ning binary black holes (BBHs) with mass-ratio 1 (mass-ratio 
4) of observed total mass ^ 2 OOM 0 and ^ SOOM©, the 
distance reach is ^ 5 Gpc (respectively, ^ 3 Gpc), with 
the largest reach of ^ 6.5 Gpc (respectively, ^ 4 Gpc) for 
^ 4 OOM 0 (see Fig. 1). The intrinsic mass (i.e., the rest-frame 
mass) of a binary M is related to the observed mass Mobs 
by Mobs = (1 + 2 ;)M, and so the intrinsic masses detected 


^ For a network consisting of two advanced LIGO and Virgo detectors, we 
compute the distance reach as the root-mean-square distance, averaged 
over the whole sky and polarisation angle, at which the network SNR is 
equal to 12. We do not average over the inclination angle, but instead we 
use a typical value of tt/S. 


at these redshifts are signiflcantly smaller than the observed 
masses. As we see from Fig. 1, subleading modes become in¬ 
creasingly important close to coalescence and their impact on 
the SNR is relevant for BBHs of total mass > 200 M©, espe¬ 
cially for asymmetric binaries whose orbital plane is inclined 
with respect to the line-of-sight. 

The increase in the distance reach brought about by the use 
of subleading modes is greater for these latter systems as com¬ 
pared to face-on, equal-mass systems where the increase is 
negligible. When spins are included, the distance reach can 
be a factor of two larger (for near maximal BH spins aligned 
with the orbital angular momentum) or smaller (for maximal 
spins anti-aligned with the orbital angular momentum) [6-8] . 
Thus, advanced LIGO and Virgo could detect BBHs in the 
hundred solar mass range with a SNR = 12 up to 2 ; ^ 2, 
depending on the mass ratio of the system and spin. 

The above mass range falls in the domain of so-called in¬ 
termediate mass black holes (IMBHs). The formation mech¬ 
anism, evolutionary history and mass function of IMBHs are 
largely unknown, as it is very difficult to observe them and 
measure their masses in the electromagnetic window. Sev¬ 
eral mechanisms have been proposed for their birth and evo¬ 
lution [8-11]. There is now substantial evidence that galactic 
nuclei contain massive BHs of millions to billions of solar 
masses but they are believed to have been seeded by lighter 
BHs of hundreds or thousands of solar masses (for a review, 
see, e.g. Refs. [10, 12]). While there is also firm sup¬ 
port for the existence of stellar mass BH candidates [13], the 
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IMBH population seems to be missing and there is only indi¬ 
rect evidence of their existence. For example, it is suspected 
that IMBHs could be responsible for ultra-luminous X-ray 
sources. While not all such sources are believed to host an 
IMBH [14], some of them do show evidence of BHs of tens 
to hundreds of solar masses. These include a stellar mass BH 
of < ISMq in NGC7793 [15], a more massive 2OM0-3OM0 
BH in MlOl ULX-1 [16] and a - 4OOM0 IMBH in M82 [17]. 

At present we do not know of any IMBH binaries 
(IMBHB). However, astrophysical scenarios of their forma¬ 
tion have been proposed in the literature, which include hi¬ 
erarchical growth of black holes at galactic nuclei by accre¬ 
tion of gas, stars and compact objects (i.e. neutron stars and 
black holes) and dynamical capture of smaller black holes by 
nuclear black holes in stellar clusters. Hierarchical models 
of structure formation predict that supermassive BHs found 
in galactic nuclei might initially be IMBHs that grow to 
their current size by accreting gas and merging with other 
IMBHs [9, 10, 12,18-20]. In such a scenario we might expect 
mergers of IMBHBs when the Universe began assembling the 
large structure at high redshift (z ^ 10-20). Such mergers 
might have continued in the local Universe, but it is very dif¬ 
ficult to compute merger rates as we do not fully understand 
the initial conditions for IMBHs (mass function of seed BHs 
and their spins), their binaries (orbital parameters at forma¬ 
tion and population as a function of mass ratio), or the process 
by which they grow (accretion of gas and merger with other 
BHs). 

Besides growing their mass by dynamical capture in stellar 
clusters, massive BHs may form from the collapse of massive 
stars and until recently both observations and theoretical argu¬ 
ments suggested that stars above 15OM0 do not form at non¬ 
zero metallicity. However, recent observations of several stars 
with current masses larger than I 5 OM 0 in the R136 region 
of the Large Magellanic Cloud triggered a re-analysis [8] of 
the possibility that very massive BHs can have stellar origin. 
Ref. [8] found that very massive stellar-origin BHs with mass 
larger than IOOM 0 can form only in low-metallicity environ¬ 
ments (i.e., X < 0.1-0. 4 X 0 ), if the initial mass function ex¬ 
tends above 5 OOM 0 and pair-instability supernovae do not de¬ 
stroy stars with mass above 5OOM0. Moreover, the formation 
of close massive BH binaries requires that the very massive 
stars above 5 OOM 0 expand by a factor of 2 and go through 
and survive a common envelope phase. If these requirements 
are met, then massive BH binaries are expected to have mass 
ratios of at most a few, spins primarily aligned with the orbital 
angular momentum, and negligible eccentricity when they en¬ 
ter the advanced LIGO band [21]. If the above requirements 
are not met, then they will have too wide a separation to co¬ 
alesce within a Hubble time. However, other phenomena in 
dense stellar environments (e.g., cluster binary-single interac¬ 
tions) and in low-density field populations (e.g, Kozai mech¬ 
anism in triple systems) can lower the coalescence time of 
wide massive BH binaries. The investigation carried out in 
Ref. [8] concluded that on the order of a few massive BH bi¬ 
naries of stellar-origin could be observed by advanced LIGO 
and Virgo. However, due to astrophysical and theoretical un¬ 
certainties, the number of detections per year can be as high 



FIG. 1. We show the distance reach as a function of the observed 
total mass Mobs for several values of the binary mass ratio q. The 
reach is computed using a detector network consisting of two ad¬ 
vanced LIGO interferometers and advanced Virgo, using a network 
SNR =12. Detector sensitivities are given by the advanced LIGO 
and Virgo design curves [22] ; the advanced LIGO design is the zero- 
detuned high-power (ZDHP) noise curve. The right y-axis shows 
the redshift computed assuming cosmological parameters measured 
by the Planck satellite [23]. The continuous curves use inspiral- 
merger-ringdown waveforms with the most dominant five modes 
(EOBNRv2HM), while the dashed curve only includes the (2, 2) mode 
(EOBNRv2). The values here are numerically averaged over sky 
location and polarization, with fixed orbital phase and inclination: 
{0, ^jn} = {0,7r/3} rad. The coalescence time is also fixed to a 
GPS time of tc = 1000000008 s, corresponding to Sept. 14, 2011 
01:46:33 UTC. 


as hundreds or as low as zero. 

In this paper we use state-of-the-art waveform models to 
explore how well advanced GW detectors can measure the 
physical parameters of an IMBHB. Signals from IMBHB co¬ 
alescences have several important features that should be in¬ 
corporated in a study of how GW observations will help to 
measure the parameters of such systems. First, as several pre¬ 
vious studies have already pointed out (see, e.g.. Ref. [24] 
and references therin), in advanced GW detectors, the plunge, 
merger and quasi-normal-mode ringdown phases of evolution 
contribute significantly to the detectors’ distance reach if the 
binary has a total mass larger than about ^ 5 OM 0 . This means 
that we must use the full signal, that is not only the adia¬ 
batic inspiral phase, but also the merger and ringdown por¬ 
tions. Second, binaries formed in the field will most likely 
have negligible eccentricity [25] as they enter the sensitivity 
band of advanced detectors and can be assumed to trace quasi¬ 
circular orbits. For binaries undergoing dynamical capture or 
Kozai mechanism in star clusters, advanced LIGO and Virgo 
might detect mild eccentricities [26], if Mobs ^ 1 O- 2 OM 0 . 
For massive BHs, we expect negligible eccentricities when the 
binary enters the detector band. Indeed, for a fixed mass ratio 
and speed at infinity, the pericenter distance at capture is pro¬ 
portional to the total mass [21]. Thus, the frequency at capture 
is inversely proportional to the total mass. As a result, larger 
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total masses result in lower capture frequencies and thus cir¬ 
cularize more by the time the binary gets to a fixed frequency, 
such as 10 Hz. 

In this study we assume our systems to have zero eccentric¬ 
ity. Thus the gravitational wave emission in comparable mass 
binaries will be dominated by the = 2,m = 2) mode at 
twice the orbital frequency, at least until merger. Asymmetric 
systems with unequal masses, nevertheless, emit radiation at 
other multiples of the orbital frequency or subleading modes 
(see, e.g.. Sec. 10.4 in Ref. [27]). As shown by several au¬ 
thors [28-35] these subdominant modes can be important in 
the inspiral phase in improving the accuracy with which pa¬ 
rameters are deduced from GW observations, especially when 
the mass ratio of the binary is large. The amplitude of those 
subleading modes grow more and more toward merger [5, 36- 
38]. As a consequence, relevant properties of the progeni¬ 
tor binary can be recovered from the relative amplitudes of 
the subdominant modes excited in the BH remnant [39-44] 
and tests of general relativity [45, 46] can be carried out 
when those subdominat modes are included during merger 
and ringdown. Moreover, for the purposes of signal candidate 
detection in template bank searches. Ref. [47] showed that 
when constructing banks for BBH searches, the inclusion of 
subdominant modes yielded improved sensitivity for systems 
with Mobs ^ IOOMq and ^ ^ 4 (where q = mil m 2 > 1 is 
the mass ratio). Third, waveforms from nonspinning BH bi¬ 
naries on quasi-circular orbits are very simple chirp-like sig¬ 
nals, with monotonically increasing frequency and amplitude. 
However, BH spins can cause amplitude and phase modula¬ 
tions, so one must ideally include spin effects in the waveform 
model, unless the IMBH formation scenario strongly suggests 
nonprecessing or negligible spins [8], if the BHs grow their 
mass through multiple mergers in stellar clusters. In this pa¬ 
per, however, we will limit ourselves to nonspinning BH bina¬ 
ries as waveforms that include both spin effects and sublead¬ 
ing modes are not yet available. 

The rest of the paper is organized as follows. In Sec. II 
we describe the inspiral, merger and ringdown template fam¬ 
ily used in our analysis, its parameters and the main features 
introduced by the subdominant modes. In Sec. Ill we review 
the basics of Bayesian inference, the sampling technique that 
we use (i.e., nested sampling) and the priors employed in our 
study. In Sec. IV we discuss how the Bayesian evidence of 
the GW signal can be used to confirm detection, how the false 
alarm probability can be obtained from the Bayesian evidence, 
and how the Bayesian evidence changes depending on the in¬ 
clusion of the subdominant modes. In Sec. V we discuss how 
Bayesian parameter measurement depends on the binary’s to¬ 
tal mass, mass ratio, inclusion of subleading modes and pri¬ 
ors, and compare our study to previous ones. We also discuss 
the astrophysical implications of these measurements of IMB- 
HBs. Finally, in Sec. VI we draw our main conclusions. 


II. WAVEFORMS 

In this section we will discuss the waveform family used in 
this study and the parameters used to describe the signal as ob¬ 


served by a detector. In particular, we discuss the importance 
of the subleading modes and the merger and ringdown phases 
of the signal for IMBHBs. We demonstrate this by first plot¬ 
ting the signal as observed by an advanced detector and dis¬ 
cuss how the SNR is accumulated as a function of time. We 
will also plot the signal power spectrum and highlight the rel¬ 
evance of subleading modes for unequal-mass systems whose 
orbital plane is inclined with respect to the line-of-sight. 


A. Inspiral-merger-ringdown waveforms and parameters 

In this study we employ nonspinning waveforms con¬ 
structed within the effective-one-body (FOB) formalism [48, 
49] and calibrated to highly-accurate numerical relativity 
(NR) simulations having typical length of 30-40 GW cycles 
and mass ratios q < 6 [5]. More specifically, we use the 
EOBNRv2HM code in the LIGO Algorithm Library (LAL) 
[50] to generate the FOB waveform model in Ref. [5], which 
includes four subleading modes, namely the (/,m) = (2,1), 
(3,3), (4,4) and (5, 5) modes, as well as the leading (/, m) = 
(2,2) mode [These modes come from the decomposition of 
the GW signal h = — thx into —2 spin-weighted spheri¬ 

cal harmonics - 2 ^/^ [5].] Consequently, GW signals that we 
study contain the first five harmonics of the orbital frequency. 
During most of the early inspiral phase only the (2,2) mode, 
and, in particular, its Newtonian amplitude, will be the domi¬ 
nant component. Other harmonics and PN corrections become 
increasingly important as we get close to merger. The effect 
of these higher modes is especially relevant when the binary 
in question has its merger and ringdown frequencies in the 
most sensitive part of a detector’s response. The ringdown 
frequency of the final remnant of binaries consisting of non¬ 
spinning BHs of total mass 50 Mq to 500 Mq varies over the 
range 40 Hz to 400 Hz—the frequency range where LIGO and 
Virgo have the best sensitivity, and this provides the motiva¬ 
tion for our choice of masses used in this study. 

A Markov-chain Monte-Carlo (MCMC) study [32] demon¬ 
strated that the BOB waveforms of Ref. [5] are indistinguish¬ 
able from the NR waveforms [52] used to calibrate them up to 
SNR = 50 for advanced LIGO detectors. Subsequent investi¬ 
gations carried out in Ref. [53] verified the accuracy of these 
nonspinning EOB waveforms in the entire sensitivity band of 
advanced LIGO detectors and suggested that the EOB model 
be accurate even outside the region of calibration, i.e. when 
g > 6. This expectation was verified by the very good agree¬ 
ment found against the g = 10 NR waveform of 20 GW cycles 
in Ref. [54] and, especially, against the g = 7 NR waveform 
of 350 GW cycles recently produced by the SXS collabora¬ 
tion [55]. 

Since in our study we consider BBHs in quasi-circular or¬ 
bits with negligible spins, the system can be described by nine 


^ An EOB model with different parametrization was subsequently calibrated 
to the same set of NR waveforms used in Ref. [5] and provides two sub¬ 
leading modes (/, m) = (2,1) and (3, 3) [51], besides the leading one. 
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FIG. 2. We display EOBNR waveforms with subleading modes used in this study. We plot the plus polarization of each waveform: original 
(dotted blue), and normalized by the PSD (and then scaled to have equal amplitude at t = 0 as original) (solid red). The plots in the top row 
have Mobs = SOM© and in the bottom row have Mobs = dOOM©; the left column has g = 1 and the right column has g = 4. The time axis 
is scaled by the observed total mass of the system. Vertical lines are at intervals of 10% of signal power with the right most line at 99%. All 
systems are observed at an inclination of ^jn = tt/S rad and a distance of = 1 Gpc. 


parameters: 6 = {mi, m 2 , ^jn, 0}- The pa¬ 

rameters mi and m 2 are the masses of the individual BHs. 
From these quantities we define the total intrinsic mass, M = 
(mi +m 2 ), the total observed mass, Mobs = M (1-h^) with z 
being the redshift, the mass ratio, q = mi/m 2 > 1, the sym¬ 
metric mass ratio, r] = mim 2 /(mi + m 2 )^ = ^/(l + q)‘^, 
the intrinsic chirp mass, M. = and the observed chirp 

mass, Alobs = M The parameter dL is the luminosity 

distance and when combined with the declination 5 and right 
ascension a, it defines the sky location of the binary. The time 
of the peak in the (2, 2) mode of the waveform, as measured at 
the geocenter, is given by tc\ this serves as an approximation 
of the merger time. The angle 6 >jn measures the inclination 
of the binary’s total angular momentum J (equal to the orbital 
angular momentum, L, as the holes are non-spinning) with 
respect to the line of sight from the detectors n (geocenter). 
The polarization 7 /) and the phase (j) provide the additional Eu¬ 
ler angles necessary to describe the rotation from n to J. 


B. Accumulation of SNR in waveforms 

In this paper we are interested in studying the effect on pa¬ 
rameter estimation of binary systems with larger and larger 
total masses. When we maintain a constant SNR and increase 


the binary’s total mass, the signal moves downward in fre¬ 
quency space, resulting in more power from the merger and 
ringdown portions as opposed to the inspiral. The merger oc¬ 
curs at approximately the frequency of the last stable orbit 
(LSO), which in the case of a Schwarzschild BH is /lso — 
4400 (M©/M)Hz. We begin our integration at /min = 10 Hz, 
so for systems with Mobs ^ dOOM© there will be little power 
from the inspiral portion of the waveform. As the inspiral evo¬ 
lution is dominated by Alobs and the merger and ringdown are 
dominated by Mobs, we expect the character of the parameter 
estimation to transition from one to the other (see Sec. V). 

To demonstrate this expectation, we compare in Fig. 2 two 
waveforms with Mobs = 50M© and 400M©. We show 
both the original waveform and the waveform normalized by 
the advanced LIGO zero-detuned high-power (ZDHP) power 
spectral density (PSD), S'n(/); the normalized waveform has 
been re-scaled so that it has the same amplitude at f = 0 
as the original. The comparisons clearly show that for the 
(Mobs, q) = (50M©, 1) waveform, 90% of the power (SNR^) 
has been accumulated during the inspiral and that the merger 
and ringdown play only a minor role. On the other hand, 
for the (Mobs, q) = (400M©, 1) waveform, only 10% of the 
power is collected during the inspiral and now the merger and 
ringdown are predominant features. In Fig. 3 we show the 
amplitudes of the waveforms in frequency space in compari- 
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FIG. 3. We show the Fourier-domain amplitudes of the waveforms 
displayed in Fig. 2. Extra structure from subdominant modes can 
be clearly seen in the non-equal mass cases. For comparison, we 
also display the advanced L IGO ZD HP and advanced Virgo design 
amplitude spectral density \/Sn{f)- For all waveforms we use dL = 
1 Gpc and ^jn = tt/S rad. 


son to the advanced LIGO ZDHP and advanced Virgo design 
\/Sn{f)- We can see in Fig. 3 that for the higher-mass wave¬ 
forms, the entire inspiral signal is strongly down-weighted by 
the rising of the amplitude spectral density with the merger 
occurring as \/ Sn{f) reaches a minimum. For the lower- 
mass waveforms, much of the inspiral is in the frequency band 
where the amplitude spectral density is at or near minimum, 
thereby allowing this part of the waveform to dominate. 


III. BAYESIAN INFERENCE 


This section will provide a short background to Bayesian 
inference. We will focus on the application of Bayesian anal¬ 
ysis to the problems of detection, parameter estimation and 
model selection that will be used in Secs. IV and V in the 
context of GW observations of IMBHBs. After a brief intro¬ 
duction to the basics of Bayesian methods, we will discuss a 
specific technique called nested sampling that is used to effi¬ 
ciently compute Bayesian evidence, followed by a description 
of our choice of prior probabilities for various parameters and 
how we compute the likelihood function. 


A. Basics of Bayesian methods 


where Vy{9\V^1-L) is the posterior probability distribution of 
parameters used for making inferences about which signal pa¬ 
rameters 6 best fit the data and what the corresponding cred¬ 
ible regions are; Pr(D|0,H) is the likelihood of obtaining 
the data given the specific model and parameters, for which 
we use the shorthand C{9)\ Vi{9\l-L) is the prior probability 
of the parameters for the model that represents our knowl¬ 
edge of these values before looking at the data (a priori)', and 
Vr{V\l-L) is the Bayesian evidence, which is commonly abbre¬ 
viated as Z. The evidence is the factor needed to normalize 
the posterior distribution and can therefore also be expressed 
as 


Z= [ C{d)Vr{d\'H)d^d, (2) 

JQ 

where N is the dimensionality of the parameter space. Since 
Z is independent of the parameters, it can be safely ignored 
for parameter estimation problems, but it is still useful in 
model comparison. 

When comparing two models, Ho and Hi, one can write 
their relative probabilities as 

Pv{n,\v) _z^ Pvin,) 

PT{no\V) ZoPrino)' 

where we used Bayes’ theorem again and have cancelled out 
Pr(D) and substituted in Zi = Pr(D|Hi) as appropriate. The 
relative probability of the two models is thus the ratio of their 
Bayesian evidences multiplied by the relative probability prior 
to considering the data. In our analysis, we will take the lat¬ 
ter to be 1 and consider the ratio of the Bayesian evidences, 
which is called the odds ratio. In the problem of signal de¬ 
tection, Ho can be considered the noise-only model while Hi 
is the signal-plus-noise model. Therefore, an odds ratio much 
greater than 1 indicates a strong belief in the presence of a sig¬ 
nal. This method naturally incorporates Occam’s razor, such 
that more complicated models are penalized and must suffi¬ 
ciently improve the fit to the data to be favored. 


B. Nested sampling and MultiNest 

Nested sampling [56] is a Bayesian inference technique de¬ 
veloped for the calculation of the evidence, through which 
posterior probability samples are produced as a by-product. 
This is done by transforming the A^-dimensional integral for 
Z into a 1-dimensional integral over the prior volume. We 
define the prior volume X by dX = FT{0\'H)d^6. We can 
therefore write the prior probability volume enclosed within a 
contour (in parameter space) of constant likelihood A as 


Bayesian inference provides a statistically rigorous method 
of measuring the probability distribution of a set of parameters 
6 given a model or hypothesis H and a set of data V. Bayes’ 
theorem states that The evidence integral of Eq. (2) can be re-written as 


X{X) = ( Pi:{e\'H)d^e. 

Jc(0)>\ 


Pr:{e\V,n) 


Pv{v\e,H) Pv{e\n) 

Pr(P|-H) 


( 1 ) 


^ = 


f 


L{X)dX, 


(4) 


(5) 
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where L{X) is the inverse of Eq. (4) (returns the likelihood 
at which a prior volume of X is enclosed) and is a mono- 
tonically decreasing function of X (i.e. more prior volume 
implies lower likelihood contour bound). If we can evaluate 
likelihood values Li = L{Xi) such that Xi is a sequence of 
monotonically decreasing values, the evidence can be com¬ 
puted as a simple sum 


M 

^ = 'L 

i=l 

Here, the Wi are weights which can be taken from a simple 
trapezium rule such that Wi = ^{Xi_i — X^+i). 

The individual prior weight of each sampled point can also 
be estimated from the sequence of Xi values. This may be 
combined with the computed likelihood for that point and the 
evidence to produce a final posterior probability for the point. 
The full sequence of points can then be re-sampled accord¬ 
ingly to the points’ individual probabilities to produce a set of 
samples from the posterior. 

Nested sampling operates by starting with an initial set of 
dive’ points sampled from the prior distribution. Iterations 
are then performed whereby the point with lowest likelihood 
value is removed from the live point set and a new point is 
sampled from the prior with the restriction that it has higher 
likelihood than the point just removed. This removal and re¬ 
placement is continued until a stopping condition is reached 
(e.g., a tolerance on the evidence calculation). The difficult 
task here lies in the efficient sampling of new points under 
this restriction. As the likelihood contour moves upwards, the 
volume of the prior within that contour will decrease to very 
small values, making direct sampling of the prior very inef¬ 
ficient. The MultiNest algorithm [57-59] addresses this 
by enclosing the live points in clusters of ellipsoids. A new 
sample can then be made from the ellipsoids very quickly and 
as they shrink along with the live points, they create effec¬ 
tive likelihood contours to be sampled from, thereby greatly 
increasing the sampling efficiency. The ellipsoids can be dis¬ 
tributed to enclose degenerate and multimodal distributions, 
making this approach very robust. 

MultiNest is implemented within BAMBI [60], which 
is linked with the LALInference [61] code of LAL. The 
lalinference_b amb i sampler is used for analysis of sim¬ 
ulated signals in this study. 


C. Priors used 

The priors used in this analysis are fiat in the component 
masses, with {mi, m 2 } G [10, 600] Mq and mi > m 2 . In 
Section VD, we assess the effect of changing the prior by 
implementing an alternative mass prior that is fiat in the log 
of chirp mass for M G [2.45,435.275] Mq and fiat over r] G 
[0.03,0.25]. 

In both setups, the source location prior is uniform in vol¬ 
ume, thus proportional to for di, G [100 Mpc, 10 Gpc] and 
fiat in sin(J) and a for G [0, tt] rad and a G [0, 27r) rad. We 
use a prior fiat in coalescence time that is centered on the true 


value with Ate G [—0.1, 0.1] s. The orientation angles are as¬ 
sumed to be isotropically distributed, thus fiat in 'ip, 

and (j) for G [0, tt] rad, ^ G [0, tt) rad, and (j) G [0, 2tt) rad. 


D. Likelihood function 


In general, the data obtained from adavanced LIGO and 
Virgo detectors is the sum of signal, h, and noise, n, 

d = h -f n. (7) 


The signal in a given detector is given by 

h = F+(q:,( 5, V’)h++ Fx(a,(5, V’)hx, (8) 


where h+^x are the two independent GW polarizations and 
X are the antenna response functions [62] that de¬ 

pend on the source location and polarization. The antenna 
response is slowly varying in time due to the rotation of the 
Earth, but this effect is small for the short duration of the sig¬ 
nals considered in this study (< 2 min). 

The noise is modeled as independent and Gaussian in each 
frequency with a mean of zero and variance given by the de¬ 
tector’s PSD. Therefore, the probability of a data stream, d^, 
in detector a containing a given signal, h(0), is given by the 
probability of the resulting noise realization, n = d — h(0). 
This is given by the product (sum in log-space) of the proba¬ 
bility of the noise for each frequency bin [61]: 


\ogCa{e)=\ogY>l{da\Hs^e,SM)) = 

b , , (TTTSMi)\ 

2 A SMi) ® V 2 J 


(9) 


where T is the segment length, the tilde indicates the discrete 
Eourier transform of the function, and i is an index over fre¬ 
quency bins. The noise power spectral density Sn{f) will vary 
from detector to detector and here we use the ones at design 
sensitivity for advanced LIGO and Virgo [22]^. Hs indicates 
that we are using the signal model that assumes a signal is 
present; this will be compared to the noise-only model, Hn, 
where h = 0. The final likelihood is the product of likeli¬ 
hoods from the individual detectors. 


log£(0)=^log£aW. (10) 

a 


To simulate the sensitivity for advanced LIGO and Virgo de¬ 
tectors, we use a minimum frequency of /min = 10 Hz. In or¬ 
der to include the highest ringdown mode for the lowest pos¬ 
sible total mass system, we use a sampling rate of 4096 Hz, 
giving a Nyquist frequency of /Nyq = 2048 Hz for the upper 
bound of our likelihood sum. A segment length of 128 s en¬ 
sures that no waveforms are cutoff in-band. 


^ See https://dec.ligo.org/LIGO-P12 0008 7-vl9/public 

for PSD data files 
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FIG. 4. False alarm probabilities and rates computed using Bayesian 
inference for signals buried in Gaussian noise as a function of the 
injected signal’s SNR. The solid line is the median FAP/FAR and the 
shaded area covers the range from minimum to maximum values of 
FAP/FAR. Confident 5-cr significant detections can be calimed for 
SNR > 9 but SNRs ~ 12 are needed in reality for the same signifi¬ 
cance in real data that is often non-Gaussian and nonstationary. 


IV. BAYESIAN DETECTION 


It is computationally infeasible to perform a Bayesian anal¬ 
ysis over the entire detector data set, for all signal types, at 
all times. Therefore, alternative analysis pipelines are used 
to first produce candidate triggers for follow-up analyses us¬ 
ing Bayesian inference over a small data set and parameter 
space [63, 64]. In the search for BBHs, a discrete bank of 
template waveforms [65-67] is used to perform matched filter 
analysis of the data. The matches that cross a pre-set thresh¬ 
old are ranked by a re-weighted SNR [63, 64] and their sig¬ 
nificance is measured by comparison to the estimated back¬ 
ground (i.e. noise-generated) triggers. This is a frequentist 
method of detection (and significance measurement) and is 
very useful for generating triggers which then receive a more 
detailed follow-up with Bayesian analysis and other tools. It 
is the first step in identifying and confirming a GW signal with 
LIGO and Virgo. 

In Bayesian inference, one can make claims on the presence 
of a signal in data by means of model comparison. This is not 
the same as trigger/candidate finding, but rather looking at the 
evidence that the candidate is indeed a real GW signal. We 
compare the signal-plus-noise and noise-only models, Hg and 
Hn, and their respective evidences, Zg and Z^. The probabil¬ 
ity that random noise would produce an evidence ratio Zg /Zn 
is the false alarm probability (FAP). This is given by 


FAP = 


1 

1 + 


( 11 ) 


The relative log-evidence (iog Zg — log Zn) is output by 
LALInference. The FAP can be converted into a false 
alarm rate (FAR) by dividing it by the length of the time win¬ 
dow [68]. This accounts for the amount of time in which we 
searched for a signal in the data and assumes that all such time 


intervals are independent. 


FAR 


1 + 


( 12 ) 


The time window used in our prior is At = 0.2s. 

We compute the Bayesian evidences of the signal plus 
Gaussian noise model for all simulated signals using in¬ 
jections at SNRs ranging from 6 to 18. Total observed 
masses ranged from SOMq to SOOMq, mass ratios were q = 
{1.25,4}, and inclinations were 6 >jn = {0, tt/S, 7r/2}. In 
Fig. 4 we show the median FAP and FAR calculated over 
all signals as a function of the SNR. The shaded area cov¬ 
ers the range from minimum to maximum computed FAPs. 
Confident detection can be claimed for a network SNR > 9, 
as this corresponds to a detection outside of the ±6-cr re¬ 
gion (FAP < 10“^). This should be taken with a grain of 
salt, however, as real data will contain non-Gaussian and non¬ 
stationary noise features that will need to be addressed with 
noise modeling [69, 70]; SNRs ^ 12 are needed in reality for 
the same significance in real data that is often non-Gaussian 
and nonstationary. 

In systems where the subleading modes contribute signif¬ 
icantly to the SNR, not including them can result in not re¬ 
covering the full power of the signal. This means that using a 
(2, 2) -only template (EOBNRv2) will yield a lower Bayesian 
evidence than a more complete template (EOBNRv2HM) for 
the same signal. This loss in evidence will lead to a greatly in¬ 
creased FAP and FAR. A model comparison between the two 
will favor the complete model when these modes are signifi¬ 
cant - inclined systems with larger mass ratios. Fig. 5 shows 
that in these cases the EOBNRv2HM waveform model includ¬ 
ing subleading modes will be strongly favored; when sublead¬ 
ing modes contribute little SNR, neither waveform model is 
strongly favored over the other (slight preference for (2,2)- 
only when 6 >jn = 0 and slight preference for subleading 
modes when 6 >jn > 0 and q is close to 1). 

The importance of including subleading modes in search 
pipelines was investigated in Ref. [47]. They found that us¬ 
ing waveforms with subdominant modes increases the sensi¬ 
tive region only for high total masses (Mobs ^ IOOMq) and 
asymmetric (q > 4) IMBHBs. Furthermore, they found that 
the most significant gains are in regions of the parameter space 
with the lowest expected event rates. Although the study of 
Ref. [47] was limited to component masses rrii < 2 OOM 0 
and Mobs < SGOMq, we can expect the trends to continue 
for larger masses. The result that subleading modes are sig¬ 
nificant in detection only for asymmetric and large total mass 
systems is consistent with our findings described in this sec¬ 
tion. 


V. MEASUREMENT 


After the detection of a GW signal from a binary system, we 
perform parameter estimation analysis, which involves pro¬ 
ducing a sufficient number of samples from the posterior dis¬ 
tribution so that we are able to measure peaks and analyze 
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FIG. 5. We show the difference in log-evidence between using 
EOBNRV2HM and EOBNRv2 templates for recovering E0BNRv2HM 
signals. In cases with significant contribution from subleading modes 
the EOBNRV2HM model is strongly favored. 

correlations and degeneracies. Since we expect our first detec¬ 
tions to be at just above threshold, all analyses in this section 
- unless otherwise stated - use injected signals with a network 
SNR of 12. This is achieved by adjusting the distance of the 
signal to obtain this exact value. 

In the following sections we discuss our ability to perform 
parameter estimation under varying conditions. We estimate 
the statistical uncertainty and bias in the measurement of sig¬ 
nal parameters; these are the width of the posterior distribu¬ 
tion and the distance between the peak and the true values, 
respectively. In creating the data to be analyzed, no noise re¬ 
alization is added. This eliminates additional uncertainty and 
bias introduced by a random noise realization; zero noise is 
the most probable realization. This is different from averag¬ 
ing over many noise realizations, as the latter would result in 
increased uncertainty even as the biases cancel out (and would 
also require many more runs to be performed). 

Results presented are predominantly foYq = {1.25,4} sys¬ 
tems. Analyses were also performed where injected wave¬ 
forms had g = {1,2, 3}; we found the results to be consistent 
with those discussed here. We limited < 4 for injected sig¬ 
nals due to the increased computational cost for higher mass- 
ratio waveforms. 


A. Measuring variance with increasing binary’s total mass 

In our first set of comparisons, we study the effect of the 
total mass of the system on the estimation of the source’s in¬ 
trinsic parameters. Specifically, we investigate the statistical 
errors on the measurement of Mobs. M,r], mi, and m 2 . Pos¬ 
terior distributions over (Mobs,Mobs)-space and (mi,m 2 )- 
space are shown for various values of Mobs (always using 
q = A and 6 >jn = tt/S) in Fig. 6 . These are summarized in 


the left panel of Fig. 7, which shows the relative widths of the 
95% (±2-cr) credible intervals (i.e., (^ 97 . 5 % — X 2 .^%)/Xi^ue) 
for the various mass parameters. At the lower mass end 
(Mobs = 5 OM 0 ), uncertainty is low due to how well the chirp 
mass Mobs is measured from the inspiral phase of a wave¬ 
form. As the total mass increases, the uncertainty increases 
and for Mobs ^ 15OM0 the uncertainty in Mobs is similar to 
or greater than that in Mobs- This change is due to less inspi¬ 
ral signal being present in the most sensitive band of the de¬ 
tector; the ringdown is predominantly dependent on Mobs and 
therefore this parameter is measured more accurately. How¬ 
ever, the inspiral measures Mobs better than the ringdown 
measures Mobs, so the resulting uncertainty is larger. Above 
Mobs = 3 OOM 0 , the uncertainty decreases slightly; this is 
due to the ringdown matching up better with the minimum 
of the advanced LIGOWirgo PSD and the subleading modes 
moving into more sensitive regions of the PSD. 

When these same systems are face-on (6 >jn = 0) or have 
lower q (more equal component masses), the uncertainties for 
the (2, 2)-only waveform closely resemble those for the wave¬ 
form with all modes, just a little larger. This can be seen in 
the right panel of Fig. 7 and it is what would be expected for 
systems with little contribution from the subleading modes. In 
all cases, the templates that include subleading modes of radi¬ 
ation have lower uncertainty than those using only the (2,2) 
mode. Additionally, for inclined and asymmetric systems, 
as Mobs increases from 3 OOM 0 to 5 OOM 0 , the uncertainty 
when using (2, 2)-only templates grows significantly while 
that from using templates with subleading modes slightly de¬ 
creases. This is due to the fact that the subleading modes pro¬ 
vide information about the mass ratio of the system in their 
relative amplitudes and phases. This information contained in 
the subdominant modes breaks the model degeneracies and al¬ 
lows us to better infer the component masses as the ringdown 
phase of the waveform enters the most sensitive region of the 
PSD. These results can also be seen in Fig. 6. Being able 
to accurately measure the component masses is important in 
allowing us to make inferences on the source population of 
these massive BHs [8]. 

Thus, we observe that at an SNR of 12 uncertainties for 
Mobs can reach ^ 20-25% in asymmetric binaries while un¬ 
certainty in Mobs reaches up to ^ 50-60% (these numbers 
are ^ 17% versus ^ 22 % in more symmetric mass ratio bina¬ 
ries). 

In Appendix A, we provide summary tables of relative 95% 
credible intervals for measurements of the masses, luminosity 
distance, and coalescence time. These are given over a range 
of SNRs for two mass ratios (q = {1.25,4}) and several ob¬ 
served total masses at the inclinations of 6 >jn = { 7 r/ 3 , 0}. 


B. Measuring degeneracies with increasing binary’s total mass 

As the IMBHB systems increase in total mass, there are dis¬ 
tinct changes in the two-dimensional posterior probability dis¬ 
tributions. As the inspiral phase evolution strongly constrains 
Mobs and 77 , in that order, lower mass systems will have de¬ 
generacies that follow contours in these parameters. With in- 
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FIG. 6 . Posterior distributions of the mass estimation. All values are presented as fractional errors, i.e., {x — Xtme)/xtme’ The left column 
displays m 2 vs. mi and the right column displays Mohs vs. Mobs- The rows are of increasing Mobs from Mobs = SOM© at the top to 
Mobs = 5 OOM 0 at the bottom. For all systems, g = 4 (77 = 0.16) and ^jn = tt/S. The star indicates the point with highest log£ and the 
contours are at 50%, 90%, and 95% credible levels (inside to outside). In the left column, the solid lines are of constant Mobs and the dashed 
lines are constant Mohsl in the right column, solid is constant mi and dashed is constant m 2 . In all cases the lines intersect the true values at 
{x,y) = ( 0 , 0 ). 


creasing total mass, however, the inspiral becomes less im¬ 
portant and the merger-ringdown part of the signal contributes 
significantly or dominantly to the SNR. This is most well de¬ 
scribed by Mobs with much weaker dependence on 77 . Thus, 
we expect there to be a change in the degeneracies present in 
the mass estimation. The inspiral dependency on A^obs can 


be seen in the PN inspiral waveforms (see Ref. [71]); PN ap- 
proximants are accurate for early inspiral when the BHs are 
sufficiently far from merger. The dependency on Mobs of 
the ringdown is similarly given by the quasi-normal mode de¬ 
composition derived in Ref. [72] and implemented in the EOB 
waveform models used in our study [5]. 
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FIG. 7. Uncertainty in the measurement of different mass parameters: Mobs, M, and 77 . Results here are for a system with (left panel) q = 4 
(?7 = 0.16) and (right panel) q = 1.25 (77 = 0.247). In both cases, ^jn = tt/S and SNR = 12. The solid lines are for templates with 
subleading modes and the dashed lines are for templates with the (2, 2) mode only. The relative width is given by (xgr. 5 % — a; 2 . 5 %)/a:true. 


This change in the optimally measured parameters as Mobs 
increases can be observed in Fig. 6. In the right column, 
we show the posterior distribution of the masses, parame¬ 
terized as Mobs and A^obs. over a range of total masses for 
an asymmetric system (q = 4, 77 = 0.16) viewed at an an¬ 
gle (^jN = 7r/3). In the top row, Mobs = SOMq and we 
can see that the principal measurement is of the chirp mass 
Afobs - posterior samples and contours lie along a line of 
near-constant Alobs- As Mobs increases, at Mobs = ISOMq 
a second principal direction of degeneracy becomes evident. 
This is due to a different combination of the mass parameters 
becoming increasingly constrained relative to the others and 
realizing a new degeneracy in the measurement. 

These observations confirm what we see in the one¬ 
dimensional posteriors in Fig. 7. The chirp mass Alobs is 
initially measured to lower fractional error than the total mass 
Mobs; as Mobs increases the uncertainty grows much faster in 
Afobs than it does in Mobs- The small decrease in uncertainty 
at the higher masses is also visible as the contours shrink 
slightly. We are now able to see in Fig. 6 that this increase 
in uncertainty is accompanied by a changing of the dominant 
degeneracy in the parameters of the waveform model. 


C. Importance of including subleading modes 

As discussed previously, in addition to the leading (2,2) 
mode, the EOBNRv2HM waveform model also includes sub¬ 
leading modes (2,1), (3,3), (4,4), and (5,5), which intro¬ 
duce additional structure to the waveform and improve faith¬ 


fulness to NR waveforms. This increased structure is impor¬ 
tant as the relative amplitudes and phasing of the additional 
modes introduce information about the source masses. In the 
ringdown phase, the additional modes further constrain the 
mass and spin of the final BH. This structure creates variation 
in waveforms as initial component masses are varied, thereby 
allowing Bayesian inference to measure the masses more ac¬ 
curately and precisely as seen in Figs. 7, 8, and 9. 

The modes’ structure also contains angular dependence on 
the inclination of the system to the detector, 6 >jn, the orienta¬ 
tion of the orbit in the plane of the sky (polarization), 7 /), and 
the orbital phase of the binary, (j). This structure is present 
for the primary (2,2) mode, but the introduction of additional 
modes breaks degeneracies in the observed waveform as these 
angles vary. 

Improvements in measurement of the masses and orienta¬ 
tion angles can all be observed in the example presented in 
Fig. 8. This figure compares one-dimensional posterior dis¬ 
tributions for the case of an injected signal with Mobs = 
5 OOM 0 , q = 4, and 6 >jn = 7r/3. The solid black line shows 
the posterior distribution for the waveform template model in¬ 
cluding subleading modes (EOBNRv2HM) and the dashed line 
is for the waveform template model including only the lead¬ 
ing (2, 2) mode (EOBNRv2). In both cases, the signal injected 
into the data contained the subleading modes - in nature, all 
modes of radiation are present. 

For parameter estimation, the presence of subleading 
modes means that as the angles 6 >jn, and (j) vary, there is 
increased variation of the waveform. This variation is more 
prominent for unequal mass binaries and binaries not ob- 
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FIG. 8. One-dimensional posterior distributions for all parameters for an injected signal with Mobs = bOOM©, g = 4, ^jn = tt/S, and 
SNR = 12. Posteriors are compared for (solid) the EOBNRv2HM waveform model which includes subleading modes of radiation, (dashed) 
the EOBNRv2 waveform model which includes only the (2, 2) mode, and (dotted) the EOBNRv2HM waveform model along with a different 
prior distribution that is flat in (log(Adobs), '^) instead of flat in (mi, m 2 ). The vertical dashed red lines indicate the true (injected) value of 
each parameter. 


served face-on or face-off (6 >jn = {0 , tt}), as the subleading 
modes will have more significant contributions to the SNR. 
The increased variation allows for more accurate measure¬ 
ment of 6 >jn and breaks degeneracies in and (j) to allow these 
two angles to be measured. As Mobs and 6 >jn are strongly 
correlated with the luminosity distance dL via the amplitude 
of the waveform, measuring the former two more accurately 
means that the latter will be measured more accurately as well. 

The measurement of the coalescence time tc is offset when 
using only the leading mode; this is likely due to slight errors 
in measuring the sky position of the source and adjustments in 
order to align the peak amplitude of the waveform at merger. 

With the additional information provided by the subleading 
modes in the EOBNRv2HM model, posterior distributions for 
all mass parameters are narrower and better centered on the 
true values. Most notably, the subleading modes and their rel¬ 
ative amplitudes differentiate better between waveforms with 
the same total mass, but different mass ratios. The system 
injected in the analysis shown in Fig. 8 has very high Mobs 
so the merger and ringdown provide the majority of the SNR. 
The observed improvement is thus due to a waveform degen¬ 
eracy in the mass of the final BH that can be broken when we 
are able to measure the final mass and spin of the BH more 
precisely. These depend strongly on the mass ratio of the ini¬ 
tial components and are further realized in the relative ampli¬ 
tudes of the subleading modes. The improved measurement of 
initial mass values from using subleading modes can also be 
seen in Fig. 9, where we compare two-dimensional posteriors 
in the masses between EOBNRv2HM and EOBNRv2 waveform 
models over a range of Mobs (q = ^ and 6 >jn = tt/S). Note 


that when Mobs = IOOMq, the posteriors are nearly identical. 
However, with increasing Mobs, the late inspiral, merger, and 
ringdown become increasingly important. The parameter es¬ 
timation bias and loss of SNR is evident when using only the 
leading mode as the posteriors do not necessarily peak at or 
strongly support the true values and the confidence intervals 
are considerably larger. 


As the observed total mass increases, it is very difficult to 
measure the mass ratio q (or rf) if the source is face-on. Indeed, 
in this case the subdominant modes are not significant; we 
find that the posteriors are identical and do not change much 
with respect to the true q. However, when ^jn = 7 r/3or7r/2, 
because of the presence of the subleading modes, it is possible 
to measure g, although the infiuence of the prior is still evident 
in tending towards smaller values. In Ref. [8], the authors 
determine that IMBHBs formed from stellar-origin massive 
BHs will likely have mass ratios q < 1.25. We find that for 
values of g > 2, in more massive (Mobs ^ SOOMq) and 
inclined systems, we will be able to say that g > 1.25 with 
certainty > 90%. 


In summary, we find that the inclusion of subleading modes 
of gravitational radiation improves the accuracy and precision 
of the estimation of the source mass parameters as well as 
some extrinsic parameters, such as distance and orientation 
angles. They are significant for asymmetric and inclined bi¬ 
naries where they contribute more to the signal’s SNR. 
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FIG. 9. Posterior distributions of the mass estimation. All values are presented as fractional errors, i.e., (x — a;true)/^true. The left column 
displays m 2 vs. mi and the right column displays Adobs vs. Mobs- The rows are of increasing Mobs from Mobs = lOOM© at the top to 
Mobs = 3 OOM 0 at the bottom. For all systems, g = 4 (77 = 0.16) and ^jn = tt/S. The asterisks indicate the point with highest log£ and 
the contours are at 50%, 90%, and 95% credible levels (inside to outside). Blue contours use EOBNRv2HM as a waveform template while red 
contours use EOBNRv2, which only includes the leading (2, 2) mode. 


D. Effect of priors 

So far, in all of the analysis runs, we have used a very large 
prior on the component masses, which was flat in (mi,m 2 ) 
space. However, one could argue for other reasonable prior 
distributions on the masses. One such alternative is to use a 
prior that is flat in log(Alobs)- The quantity log(Alobs) is 
used because Alobs is a scaling factor for the waveform am¬ 
plitude and log(Alobs) is the so-called Jeffreys prior. Addi¬ 
tionally, we employ a prior that is flat in r] for the second mass 
parameter. 

We ran multiple analyses with this second prior option, 
which is flat in (log(A4obs) , v) And that even at an SNR of 
12 the strength of the signal is sufficient to render the different 
prior distribution a minimal factor. This can be seen in Fig. 8, 


where we show the one-dimensional posteriors from a single 
analysis. More speciflcally, we display in dotted lines the ID 
posteriors of a run with EOBNRv2HM waveform model that 
uses the alternative prior, to be compared with the solid lines 
from the run with the original prior. The lines are nearly iden¬ 
tical, with differences much smaller than those from using the 
EOBNRv2 waveform; these differences from the alternative 
prior will continue to decrease as the SNR is increased. 


E. Comparison to previous parameter-estimation work with 
inspiral-merger-ringdown waveforms 

In an earlier work, Ajith and Bose [30] used inspiral- 
merger-ringdown phenomenological waveform models 
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(IMRPhenomA) to perform a study similar to ours, but 
mainly focusing on understanding how uncertainties are 
reduced when merger and ringdown phases are included. 
Their study did not include subleading modes. They exam¬ 
ined the statistical error in parameter estimation as given by 
the Fisher information matrix"^ and MCMC analyses. Our 
uncertainties calculated using EOBNRv2 templates ((2,2) 
mode only) are a factor of a few larger than those given in 
Table 1 of Ref. [30] for SNR = 10, Mobs = {100, 200}Mo, 
and T] = {0.25,0.16} (e.g.. Ref. [30] quotes relative uncer¬ 
tainties of 2.39% and 3.57% for Mobs in cases with q = I 
and Mobs = IOOMq or Mobs = 200Mq; we measure 
uncertainties of 4.01% and 9.66% in these same cases). The 
discrepancy in uncertainties is partly a result of the fact that 
in Ref. [30], the authors maximize their likelihood function 
over tc and (j) and only perform a single-detector search 
using an “effective” distance that folds in sky position and 
binary inclination effects; both of these choices have the 
effect of fixing values for parameters that we allow to vary 
in our more general analysis, thus introducing additional 
uncertainty. Furthermore, we use EOBNRv2HM waveforms 
as injections, as these most closely model NR waveforms, 
and our EOBNRv2 templates differ from IMRPhenomA 
ones, the former being more faithful to NR waveforms [54]. 
Despite those differences, the quantitative measurements are 
in general agreement and we also find agreement in the qual¬ 
itative aspects of increasing uncertainty in mass parameters 
with increasing Mobs- We report our own estimations of 
the uncertainty in Appendix A, using the full EOBNRv2HM 
waveform model. 

In Ref. [32], Littenberg et al. examined systematic and sta¬ 
tistical errors of EOBNR waveforms to assess whether those 
waveforms are indistinguishable from the NR waveforms used 
to calibrate them. The authors employed both EOBNRv2 and 
EOBNRv2HM templates to recover waveforms generated by 
NR simulations and, when using subleading modes, find sys¬ 
tematic errors to be comparable to or less than statistical errors 
for mass ratios up to g' = 6 and SNRs up to 50. We find sta¬ 
tistical errors comparable to the ones of Ref. [32] for analyses 
run in common. Our results on the importance of sublead¬ 
ing modes for unequal-mass and inclined (^jn > 0) systems 
reaffirm their findings. 

Varma et al. [37] built on Refs. [32, 47]. They used as 
targets “hybrid” waveforms constructed by attaching PN in- 
spiral/EOB waveforms to NR merger-ringdown waveforms, 
and as templates EOBNRv2 waveforms. Instead of examin¬ 
ing only a few points in parameter space, they ran many sim¬ 
ulations in order to average over relative orientation angles. 
Statistical errors were computed with the Fisher information 
matrix. Confirming previous work, they found that sublead¬ 
ing modes are more important for parameter estimation when 
Mobs > 15OM0 and g > 4 and important for detection when 


The Fisher information matrix measures covariances analytically. The 
square root of the diagonal elements of the inverse of the Fisher matrix 
gives a lower bound on the standard deviation of the posterior for the pa¬ 
rameters. In the limit of large SNR, this estimate becomes exact. 


Mobs > IOOMq and g > 6 (see Fig. 1 of Ref. [37]). In 
contrast, our paper employs the full EOBNRv2HM waveforms 
and uses Bayesian analysis, and it extends the study to higher 
Mobs- 

Bose et al. [31] focused on the importance of including 
merger and ringdown phases of the waveform. They ana¬ 
lyzed the recovery of inspiral-merger-ringdown waveforms 
with PN inspiral waveforms and found that at masses as low 
as Mobs = 5 OM 0 there are serious systematic errors and an 
increase in statistical errors due to the loss in SNR. This result 
is consistent with our analysis showing the significant amount 
of signal power present in the merger and ringdown phases of 
the waveform at Mobs = 50 Mq and above. 


E Astrophysical implications 

As mentioned in Sec. I, the merger rate of IMBHBs is cur¬ 
rently highly uncertain. Consequently, detection of a single 
event will immediately confirm the existence of these systems 
and constrain their rate. In the absence of detection the upper 
limits reached could be used to rule out some of the models. 

As shown in Fig. 1, ground-based detectors will have the 
greatest distance reach for equal-mass IMBHBs of observed 
total mass of ^ 400 Mq. The distance reach for such sys¬ 
tems will be ^ 6.5 Gpc or ^ 1. We find that the comov¬ 

ing volume averaged over all source orientations and weighted 
by the antenna pattern functions of the advanced LIGO-Virgo 
network is ^ 150 Gpc^, larger by a factor ^ 1.8 than that 
in Ref. [8]. This difference can be explained because, as op¬ 
posed to Ref. [8], we consider a detector network (which in¬ 
creases the reach) and use a different SNR (which decreases 
the reach). After five years of non-observation of IMBHBs a 
rate upper limit of 4 x 10“^^ Gpc“^yr“^ can be achieved, 
which is smaller than the rates for most formation models 
discussed in Ref. [8-10]. For a binary of same total mass, 
^ 4 OOM 0 , but mass ratio g = 4, the reach when includ¬ 
ing subleading modes is smaller by a factor of 1.5 (see Fig. 
1) and the upper limit will be larger by a factor 3.375, i.e., 
1.35 X 10“^^ Gpc“^ yi*~^- Neglecting the subleading modes 
worsens the upper limit by a factor 2.4. However, Ref. [47] 
showed that in a realistic search the improvement when in¬ 
cluding subleading modes is significant only for mass ratios 
larger than r\j 

An important question in cosmology is the mass function of 
IMBHs. Routine detection of IMBHs will help us measure the 
mass function of component BHs that form merging binaries 
and this should be a proxy for the mass function of IMBHs in 
the Universe, unless IMBHBs are formed selectively from a 
sub-population of IMBHs. The component masses of a binary 
system are strongly correlated and it is not possible to measure 
them accurately while using only the dominant mode; sub¬ 
leading modes break this degeneracy, especially in the case of 
asymmetric binaries for which the mass ratio g is large, help¬ 
ing us measure the component masses more accurately. In 
particular (see Fig. 7), the 95% credible interval in the mea¬ 
surement of the heavier companion can be 10% to 25%, while 
the lighter component is measured within 10% to 125%, de- 
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pending on the total mass of the binary. These results are far 
better than what might be possible by electromagnetic obser¬ 
vations of such binary systems. Therefore, advanced detectors 
provide the most robust way of determining the mass function 
across the range of masses from 50Mq to SOOMq . A related 
question is the mass function of IMBHBs. Referring to Fig. 7, 
the total mass is determined to within a few percent in the case 
of lighter binaries of SOMq to within 15-25% for the heaviest 
systems of 500Mq that we consider. Thus, advanced detec¬ 
tors should help determine the mass function of IMBHBs. 

VI. CONCLUSIONS 

In this paper we used state-of-the-art waveform models for 
inspiral-merger-ringdown phases of evolution to estimate un¬ 
certainties in parameters of IMBHBs with total mass Mobs = 
5OM0-5OOMq and mass ratio q = 1-4. Because for these 
systems the majority of the SNR is accumulated during the 
last stages of inspiral, merger and ringdown phases, where 
subleading modes can become comparable to the leading 
(2,2) mode, we also included in the analysis four subleading 
modes, i.e., (2,1), (3,3), (4,4) and (5,5). In particular, we 
employed the EOBNRv2 and EOBNRv2HM waveform models 
in LAL. 

Using a Bayesian analysis, we found that for unequal-mass 
systems and inclined binaries subleading modes improve the 
measurement of the mass parameters and break degenaracies 
in distance and orientation angles (see Sec. VC). As the bi¬ 
nary’s total mass increases, the merger and ringdown phases 
dominate the SNR. Since for such high-mass binaries the sig¬ 
nal resembles a burst, the measurement will extract the dom¬ 
inant frequency of oscillation of the signal, which depends 
primarily on the total mass, thus the uncertainty in total mass 
becomes smaller than the uncertainty in the chirp mass (see 
Sec. V A, Fig. 7 and Tables I and II). In contrast, for lower 
total masses, the SNR accumulates over many cycles of in¬ 
spiral and the chirp mass is better measured. The presence of 
subleading modes is less crucial for comparable-mass systems 
and face-on binaries (see Sec. V A, Fig. 7 and Table I). Inclu¬ 
sion of subleading modes allows for improved measurement 
of the mass ratio for asymmetric and inclined systems (see 
Table II). Finally, as discussed in Sec. VF, GW observations 
of IMBHBs will demonstrate the existence and shed light on 
the demographics of IMBHs, even if component masses will 
be measured only with a fractional error of (several) tens of 
percent. 

Our analysis was restricted to nonspinning IMBHs and ex¬ 
plored only part of the parameter space. These limitations 
were a consequence of the fact that EOBNR waveform models 
are expensive to generate for Bayesian analyses. Higher-mass 
binaries (Mobs ^ 4 OOM 0 ) take tens to hundreds of millisec¬ 


onds to generate, while lower-mass binaries (Mobs ^ 5 OM 0 ) 
will require up to tens of seconds; computational time will 
quickly add up as 10^-10^ waveform computations are re¬ 
quired for a complete analysis. This cost is compounded by 
the requirements both for long segments to enclose the en¬ 
tire waveform while in-band (long due to the low minimum 
frequency) and a high sampling rate in order to include the 
subleading modes in the ringdown signal. Recently, reduced- 
order models (ROM) have been developed for EOBNR wave¬ 
forms, either for spinning, nonprecessing systems [73] or non¬ 
spinning, but with subleading-mode waveforms [74]. Future 
investigations could employ these faster template families. 

While completing this work, we became aware of the study 
of Ref. [75], which includes the effect of nonprecessing spins 
using a ROM built on the SEOBNRv2 template family [73, 
76], while discarding the subdominant modes. 
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Appendix A: Summary of measurements 

In this section we present tables of the 95% credible inter¬ 
vals for five of the measured parameters: Mobs^ 7^ *^obs 5 
and tc. Results are presented for two mass ratios (q = 1.25 
and g = 4) and for two inclinations (6 >jn = Trj3 and 6 >jn = 0). 
For all parameters except tc, these values are scaled by their 
true values and converted into percentages. Uncertainty in tc 
is presented in ms. Table I presents values for g = 1.25 and 
Table II presents results for g = 4. For the two inclinations, 
results are side-by-side, with 6 >jn = 0 in parentheses. 
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SNR 

Mobs — 50Mo 

lOOM© 

150Mo 

200Mo 

300Mo 

400Mo 

500Mo 

^ Afobs / V/obs 

6 

16.810(16.658) 

6.895 (6.655) 

2.607 (2.815) 

0.612 (0.610) 

0.479 (0.515) 

0.472 (0.497) 

0.400 (0.410) 

8 

16.756 (16.464) 

0.159 (0.158) 

0.280 (0.285) 

0.288 (0.287) 

0.262 (0.251) 

0.239 (0.255) 

0.213 (0.214) 

10 

0.052 (0.049) 

0.086 (0.090) 

0.180 (0.178) 

0.202 (0.201) 

0.201 (0.200) 

0.180 (0.181) 

0.155 (0.159) 

12 

0.041 (0.040) 

0.065 (0.064) 

0.130 (0.130) 

0.166 (0.161) 

0.164 (0.169) 

0.149 (0.150) 

0.132 (0.125) 

14 

0.031 (0.033) 

0.050 (0.050) 

0.100 (0.102) 

0.137 (0.135) 

0.143 (0.139) 

0.128 (0.124) 

0.113 (0.105) 

16 

0.027 (0.028) 

0.042 (0.045) 

0.080 (0.082) 

0.116(0.118) 

0.132 (0.129) 

0.109 (0.105) 

0.090 (0.088) 

18 

0.025 (0.025) 

0.038 (0.038) 

0.070 (0.074) 

0.100 (0.104) 

0.113 (0.107) 

0.097 (0.087) 

0.083 (0.075) 

Arilri 

6 

0.874 (0.872) 

0.866 (0.863) 

0.709 (0.725) 

0.488 (0.427) 

0.475 (0.484) 

0.512(0.510) 

0.520 (0.524) 

8 

0.870 (0.871) 

0.201 (0.191) 

0.177 (0.181) 

0.182 (0.198) 

0.218 (0.227) 

0.221 (0.218) 

0.203 (0.186) 

10 

0.109 (0.104) 

0.135 (0.144) 

0.126 (0.130) 

0.132 (0.131) 

0.158 (0.146) 

0.136 (0.153) 

0.131 (0.097) 

12 

0.085 (0.082) 

0.114 (0.112) 

0.106 (0.099) 

0.113 (0.107) 

0.131 (0.129) 

0.120 (0.111) 

0.094 (0.075) 

14 

0.063 (0.067) 

0.086 (0.093) 

0.090 (0.086) 

0.089 (0.092) 

0.108 (0.100) 

0.100 (0.082) 

0.065 (0.053) 

16 

0.055 (0.057) 

0.077 (0.086) 

0.077 (0.073) 

0.080 (0.078) 

0.107 (0.096) 

0.082 (0.061) 

0.054 (0.038) 

18 

0.049 (0.050) 

0.071 (0.072) 

0.070 (0.069) 

0.072 (0.073) 

0.085 (0.074) 

0.069 (0.049) 

0.045 (0.034) 

AAdobs/TWobs 

6 

16.913 (16.613) 

6.594 (6.315) 

0.995 (1.293) 

0.605 (0.620) 

0.603 (0.602) 

0.626 (0.634) 

0.627 (0.614) 

8 

16.662 (16.350) 

0.169 (0.182) 

0.312(0.324) 

0.333 (0.338) 

0.344 (0.344) 

0.347 (0.357) 

0.304 (0.292) 

10 

0.029 (0.027) 

0.091 (0.099) 

0.199 (0.200) 

0.237 (0.239) 

0.269 (0.263) 

0.240 (0.250) 

0.213 (0.197) 

12 

0.021 (0.021) 

0.065 (0.068) 

0.149 (0.146) 

0.203 (0.195) 

0.220 (0.223) 

0.203 (0.199) 

0.170 (0.150) 

14 

0.015 (0.016) 

0.049 (0.052) 

0.116(0.118) 

0.163 (0.166) 

0.192 (0.190) 

0.178 (0.157) 

0.138 (0.123) 

16 

0.013 (0.013) 

0.043 (0.043) 

0.091 (0.095) 

0.140 (0.142) 

0.180 (0.173) 

0.146 (0.128) 

0.110(0.098) 

18 

0.011 (0.012) 

0.037 (0.036) 

0.081 (0.084) 

0.125 (0.130) 

0.149 (0.137) 

0.127 (0.105) 

0.096 (0.086) 

AdL/dL 

6 

1.819 (2.438) 

1.252(1.649) 

0.843 (1.078) 

0.629 (0.835) 

0.469 (0.655) 

0.449 (0.615) 

0.510(0.658) 

8 

2.423 (3.066) 

1.470(1.982) 

1.229 (1.569) 

0.873 (1.252) 

0.673 (0.932) 

0.687 (0.911) 

0.803 (1.090) 

10 

1.215(1.616) 

1.211 (1.559) 

1.213 (1.628) 

1.179 (1.578) 

0.973 (1.319) 

0.973 (1.375) 

1.183 (1.560) 

12 

1.140(1.513) 

1.084(1.463) 

1.104(1.526) 

1.133 (1.559) 

1.164(1.576) 

1.270(1.561) 

1.390(1.461) 

14 

1.045 (1.333) 

0.980(1.399) 

1.002(1.417) 

1.030(1.467) 

1.106(1.404) 

1.234(1.245) 

1.235 (1.009) 

16 

0.970(1.418) 

0.896 (1.365) 

0.959 (1.430) 

1.004 (1.423) 

1.036 (1.307) 

1.035 (1.010) 

0.979 (0.821) 

18 

0.978 (1.386) 

0.895 (1.373) 

0.982 (1.378) 

0.991 (1.337) 

0.976(1.118) 

0.923 (0.846) 

0.808 (0.688) 

Ate (ms) 

6 

187.8 (188.9) 

174.7 (173.3) 

48.3 (56.5) 

28.0 (26.2) 

34.8 (34.4) 

44.1 (44.4) 

46.1 (48.8) 

8 

191.3 (190.7) 

6.9 (6.8) 

11.3 (11.6) 

14.4 (14.8) 

20.9 (19.9) 

26.0 (26.2) 

30.6 (29.5) 

10 

3.5 (3.2) 

4.8 (4.8) 

8.1 (7.6) 

10.9 (10.5) 

16.3 (15.8) 

21.5 (20.2) 

24.0 (22.8) 

12 

3.0 (2.7) 

4.0 (3.9) 

6.3 (6.2) 

9.1 (8.9) 

13.8(13.1) 

17.7 (16.3) 

19.7(18.1) 

14 

2.5 (2.2) 

3.4 (3.1) 

5.1 (5.0) 

7.5 (7.3) 

11.5 (11.3) 

15.2 (13.6) 

16.4 (14.3) 

16 

2.2 (1.9) 

2.9 (2.6) 

4.3 (4.3) 

6.6 (6.6) 

10.8 (9.8) 

12.6(11.1) 

14.1 (12.3) 

18 

2.0 (1.7) 

2.6 (2.5) 

3.8 (3.8) 

5.6 (5.8) 

9.1 (8.3) 

11.3 (9.6) 

12.2 (10.5) 


TABLE L Measurement uncertainties for E0BNRv2HM signals with E0BNRv2HM templates. The values are the (relative) widths of the 95% 
credible intervals from the one-dimensional marginalized posterior distributions, scaled by the true value when indicated. For all runs, g = 1.25 
(?7 = 0.247). ^jN = tt/S and ^jn = 0 results are shown side-by-side, with the latter in parentheses. 
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SNR 

Mobs = 50Mo 

lOOM© 

150Mo 

200Mo 

300Mo 

400Mo 

500Mo 

^ Afobs / ^ohs 

6 

16.206 (17.067) 

7.659 (7.722) 

3.378 (3.896) 

1.821 (2.195) 

0.966(1.091) 

0.822 (0.930) 

0.586 (0.690) 

8 

16.749 (16.519) 

0.336 (0.452) 

0.523 (0.626) 

0.417 (0.419) 

0.362 (0.354) 

0.334 (0.320) 

0.324 (0.314) 

10 

15.992 (0.047) 

0.138 (0.136) 

0.298 (0.229) 

0.310(0.260) 

0.292 (0.258) 

0.269 (0.262) 

0.263 (0.244) 

12 

0.034 (0.033) 

0.111 (0.102) 

0.130 (0.110) 

0.233 (0.175) 

0.253 (0.212) 

0.242 (0.214) 

0.226 (0.203) 

14 

0.028 (0.025) 

0.091 (0.082) 

0.098 (0.084) 

0.176 (0.115) 

0.220 (0.170) 

0.212(0.174) 

0.182 (0.167) 

16 

0.021 (0.021) 

0.078 (0.068) 

0.084 (0.068) 

0.111 (0.081) 

0.196 (0.147) 

0.185 (0.153) 

0.172 (0.144) 

18 

0.218 (0.044) 

0.072 (0.061) 

0.075 (0.058) 

0.085 (0.065) 

0.173 (0.124) 

0.163 (0.128) 

0.142 (0.122) 

Arilri 

6 

1.347 (1.348) 

1.343 (1.342) 

1.290(1.314) 

1.140(1.233) 

1.035 (1.066) 

1.063 (1.043) 

1.021 (1.021) 

8 

1.345 (1.349) 

0.817 (0.824) 

0.785 (0.826) 

0.753 (0.771) 

0.717 (0.776) 

0.693 (0.744) 

0.720 (0.775) 

10 

1.346 (0.126) 

0.478 (0.443) 

0.714 (0.623) 

0.705 (0.650) 

0.648 (0.627) 

0.666 (0.645) 

0.641 (0.645) 

12 

0.082 (0.079) 

0.344 (0.308) 

0.455 (0.389) 

0.595 (0.490) 

0.597 (0.539) 

0.594 (0.551) 

0.590 (0.553) 

14 

0.062 (0.057) 

0.272 (0.245) 

0.348 (0.295) 

0.492 (0.354) 

0.547 (0.455) 

0.541 (0.459) 

0.497 (0.473) 

16 

0.047 (0.047) 

0.227 (0.196) 

0.288 (0.237) 

0.369 (0.279) 

0.503 (0.384) 

0.502 (0.412) 

0.468 (0.408) 

18 

0.223 (0.123) 

0.203 (0.170) 

0.255 (0.202) 

0.281 (0.226) 

0.450 (0.340) 

0.445 (0.354) 

0.415 (0.358) 

AAdobs/TWobs 

6 

20.744 (22.182) 

9.828 (10.135) 

3.308 (4.837) 

1.480 (2.081) 

1.268 (1.391) 

1.276 (1.368) 

1.131 (1.119) 

8 

21.574 (21.289) 

0.613 (0.780) 

0.917 (0.934) 

0.822 (0.793) 

0.780 (0.761) 

0.763 (0.745) 

0.748 (0.762) 

10 

20.421 (0.027) 

0.164 (0.158) 

0.640 (0.521) 

0.691 (0.595) 

0.682 (0.628) 

0.666 (0.639) 

0.644 (0.623) 

12 

0.016 (0.016) 

0.109 (0.098) 

0.263 (0.208) 

0.551 (0.420) 

0.609 (0.517) 

0.597 (0.534) 

0.579 (0.528) 

14 

0.011 (0.011) 

0.080 (0.072) 

0.150 (0.130) 

0.428 (0.274) 

0.547 (0.430) 

0.545 (0.443) 

0.476 (0.446) 

16 

0.008 (0.009) 

0.064 (0.057) 

0.114 (0.101) 

0.267 (0.194) 

0.485 (0.366) 

0.479 (0.392) 

0.450 (0.378) 

18 

0.208 (0.033) 

0.056 (0.049) 

0.098 (0.084) 

0.179 (0.151) 

0.432 (0.316) 

0.424 (0.330) 

0.379 (0.330) 

AdL/dL 

6 

2.340 (3.107) 

1.792 (2.252) 

1.503 (2.007) 

1.259 (1.658) 

1.061 (1.467) 

1.119 (1.433) 

1.252(1.537) 

8 

3.197 (4.272) 

2.847 (3.761) 

2.465 (3.254) 

2.017 (2.683) 

1.631 (2.152) 

1.693 (2.014) 

1.909 (2.189) 

10 

3.869 (1.842) 

1.226(1.453) 

2.090 (2.034) 

1.989 (2.128) 

1.717 (1.780) 

1.691 (1.534) 

1.740(1.495) 

12 

1.144 (1.388) 

0.926(1.014) 

1.019(1.048) 

1.330(1.271) 

1.338 (1.256) 

1.272(1.127) 

1.263 (1.039) 

14 

0.988 (1.147) 

0.785 (0.875) 

0.750 (0.821) 

0.995 (0.941) 

1.083 (0.940) 

1.006 (0.889) 

0.926 (0.839) 

16 

0.878 (1.036) 

0.688 (0.765) 

0.656 (0.693) 

0.747 (0.719) 

0.924 (0.806) 

0.852 (0.743) 

0.803 (0.676) 

18 

2.011 (1.183) 

0.609 (0.670) 

0.609 (0.615) 

0.619 (0.634) 

0.795 (0.668) 

0.750 (0.644) 

0.657 (0.589) 

Ate (ms) 

6 

189.9 (190.2) 

185.6 (185.2) 

142.4 (166.4) 

58.8 (100.6) 

54.6 (57.4) 

65.2 (65.3) 

71.4 (66.7) 

8 

191.5 (188.9) 

11.0(11.9) 

17.3 (19.6) 

19.2(19.3) 

26.8 (25.8) 

36.7 (34.2) 

44.6 (44.7) 

10 

188.8 (4.4) 

7.0 (6.7) 

11.1 (10.0) 

14.9(13.1) 

21.3 (19.0) 

26.1 (24.2) 

32.3 (28.1) 

12 

3.6 (3.4) 

6.0 (5.3) 

8.2 (7.0) 

11.8 (9.7) 

18.0(15.5) 

22.3 (19.3) 

24.7 (21.2) 

14 

3.1 (2.9) 

5.1 (4.5) 

6.7 (5.7) 

9.4 (7.8) 

15.0 (12.0) 

19.1 (15.6) 

19.7 (16.7) 
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2.9 (2.6) 

4.5 (3.9) 

5.8 (4.8) 

7.6 (6.1) 

13.1 (10.6) 

15.8(13.1) 

17.6 (13.6) 
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9.0 (5.1) 

4.0 (3.4) 

5.1 (4.3) 

6.5 (5.4) 

11.9 (8.7) 

13.6 (10.6) 

14.4(11.6) 


TABLE IL Measurement uncertainties for E0BNRv2HM signals with E0BNRv2HM templates. The values are the (relative) widths of the 95% 
credible intervals from the one-dimensional marginalized posterior distributions, scaled by the true value when indicated. For all runs, g = 4 
(?7 = 0.16). ^jN = tt/S and ^jn = 0 results are shown side-by-side, with the latter in parentheses. 
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